#Daniel J. Mallinson
#Reproduction Code for "Agenda Instability in Pennsylvania Politics: Lessons for Future Replication
# 11-23-2015
# 01-25-2016

rm(list=ls())


#install.packages(c("foreign", "moments", "Hmisc", "EnvStats", "lmom"))

library(foreign)
library(moments)
library(Hmisc)
library(EnvStats)
library(lmom)

###############################################################################
############################## Spending Kurtosis ##############################
###############################################################################

#Load spending data for kurtosis test
spending <- read.csv("correct_spending_data.csv")

data <- spending
data <- na.omit(data) #remove blank pct_diff cells

data <- data[which(data$pct_diff<300),]

head(data) #Look at the data

specify_decimal <- function(l,m) format(round(l,m), nsmall=m)

#############################################################################
############################# Figure 1 ######################################
#############################################################################

#setEPS()
#postscript("All_Spending_Histogram.eps")
pdf("All_Spending_Histogram.pdf")
#tiff("All_Spending_Histogram.tiff", height=6, width=6, units="in", res=600)
x <- data$pct_diff
kurtosis <- specify_decimal(kurtosis(x), 2)
kurtosis
l.moment <- round(samlmu(x), digits=2)
l.moment
l.kurtosis <- l.moment[4]
h <- hist(x, breaks=400, main="", xlab="Percent Changes in All Spending", ylab = "Frequency", ylim=c(0,100), xlim=c(-100,300), col="gray37", border="gray37")
xfit <- seq(min(x), max(x), length=100)
yfit <- dnorm(xfit,  mean=mean(x), sd=sd(x))
yfit <- yfit*diff(h$mids[1:2]*length(x))
lines(xfit, yfit, col="black", lwd=2, lty=1)
abline(v=0, col="black", lwd=2)
abline(v=median(x), col="black", lwd=2, lty=2)
kurt.lab <- bquote(bold(L-kurtosis == .(specify_decimal(l.kurtosis,2))))
med.lab <- bquote(bold(Median == .(specify_decimal(median(x),2))))
text(150, 80, labels=kurt.lab)
text(150, 75, labels=med.lab)
dev.off()

#############################################################################
############################# Figure 2 ######################################
#############################################################################

#Education
x <- data$pct_diff [subset=data$major_topic==6]
kurtosis <- specify_decimal(kurtosis(x), 2)
kurtosis
l.moment <- round(samlmu(x), digits=2)
l.moment
l.kurtosis <- l.moment[4]

# Labor, Employment, and Immigration
x <- data$pct_diff [subset=data$major_topic==5]
kurtosis <- specify_decimal(kurtosis(x), 2)
kurtosis
l.moment <- round(samlmu(x), digits=2)
l.moment
l.kurtosis <- l.moment[4]

# Banking, Finance, and Domestic Commerce
x <- data$pct_diff [subset=data$major_topic==15]
kurtosis <- specify_decimal(kurtosis(x), 2)
kurtosis
l.moment <- round(samlmu(x), digits=2)
l.moment
l.kurtosis <- l.moment[4]

#Law, Crime, and Family
x <- data$pct_diff [subset=data$major_topic==12]
kurtosis <- specify_decimal(kurtosis(x), 2)
kurtosis
l.moment <- round(samlmu(x), digits=2)
l.moment
l.kurtosis <- l.moment[4]

#Environment
x <- data$pct_diff [subset=data$major_topic==7]
kurtosis <- specify_decimal(kurtosis(x), 2)
kurtosis
l.moment <- round(samlmu(x), digits=2)
l.moment
l.kurtosis <- l.moment[4]

#State Government Operations
x <- data$pct_diff [subset=data$major_topic==20]
kurtosis <- specify_decimal(kurtosis(x), 2)
kurtosis
l.moment <- round(samlmu(x), digits=2)
l.moment
l.kurtosis <- l.moment[4]

# Energy
x <- data$pct_diff [subset=data$major_topic==8]
kurtosis <- specify_decimal(kurtosis(x), 2)
kurtosis
l.moment <- round(samlmu(x), digits=2)
l.moment
l.kurtosis <- l.moment[4]

#Health
x <- data$pct_diff [subset=data$major_topic==3]
kurtosis <- specify_decimal(kurtosis(x), 2)
kurtosis
l.moment <- round(samlmu(x), digits=2)
l.moment
l.kurtosis <- l.moment[4]

#Public Lands and Water Management
x <- data$pct_diff [subset=data$major_topic==21]
kurtosis <- specify_decimal(kurtosis(x), 2)
kurtosis
l.moment <- round(samlmu(x), digits=2)
l.moment
l.kurtosis <- l.moment[4]

#Agriculture
x <- data$pct_diff [subset=data$major_topic==4]
kurtosis <- specify_decimal(kurtosis(x), 2)
kurtosis
l.moment <- round(samlmu(x), digits=2)
l.moment
l.kurtosis <- l.moment[4]

#Community Development and Housing
x <- data$pct_diff [subset=data$major_topic==14]
kurtosis <- specify_decimal(kurtosis(x), 2)
kurtosis
l.moment <- round(samlmu(x), digits=2)
l.moment
l.kurtosis <- l.moment[4]

#Transportation
x <- data$pct_diff [subset=data$major_topic==10]
kurtosis <- specify_decimal(kurtosis(x), 2)
kurtosis
l.moment <- round(samlmu(x), digits=2)
l.moment
l.kurtosis <- l.moment[4]

#Social Welfare
x <- data$pct_diff [subset=data$major_topic==13]
kurtosis <- specify_decimal(kurtosis(x), 2)
kurtosis
l.moment <- round(samlmu(x), digits=2)
l.moment
l.kurtosis <- l.moment[4]


#################################################################
############## Citizen-Agenda Correspondence ####################
#################################################################

rm(list=ls())

specify_decimal <- function(l,m) format(round(l,m), nsmall=m)

#Load PA House hearings and most important question data for correlation analysis
congruence <- read.csv(file="correspondence_data.csv")

#Seperate into two matricies: hearing and mip proportions
#Note, agriculture2[8] and science[28] are removed from the hearings matrix and mipagriculture[9] and mipscience[29] are removed from the mip matrix because there were no recorded values for MIP

hearings <- as.matrix(congruence[c(2,4,6,10,12,14,16,18,20,22,24,26,30)])
print(colnames(hearings))
colnames(hearings) <- c("Economy Hearings", "Civil Rights Hearings", "Health Hearings","Labor Hearings", "Education Hearings", "Environment Hearings", "Energy Hearings", "Transportation Hearings", "Crime Hearings", "Welfare Hearings", "Housing Hearings", "Commerce Hearings", "State Gov't Hearings") #Change labels for correlation plot
print(colnames(hearings))

mip <- as.matrix(congruence[c(3,5,7,11,13,15,17,19,21,23,25,27,31)])
print(colnames(mip))
colnames(mip) <- c("Economy MIP", "Civil Rights MIP", "Health MIP","Labor MIP", "Education MIP", "Environment MIP", "Energy MIP", "Transportation MIP", "Crime MIP", "Welfare MIP", "Housing MIP", "Commerce MIP", "State Gov't MIP") #Change labels for correlation plot
print(colnames(mip))

#Create correlation matrix
correlation <-cor(mip,hearings)
print(correlation)

#Run correlation test to determine p-values
library(Hmisc)
rcorr(mip, hearings, type="pearson")

#############################################################################
############################# Figure 3 ######################################
#############################################################################


#Visualize correlation
library(corrplot)

cor.mtest <- function(mat) {
	mat <- as.matrix(mat)
	n <- ncol(mat)
	p.mat <- matrix(NA,n,n)
	diag(p.mat) <- 0
	for(i in 1:(n-1)){
		for(j in (i + 1):n){
			tmp <- cor.test(mat[,i], mat[,j])
			p.mat[i,j] <- p.mat[j,i] <- tmp$p.value
		}
	}
	colnames(p.mat) <- rownames(p.mat) <- colnames(mat)
	p.mat
}

p.mat <- cor.mtest(cbind(mip,hearings))
p.mat <- p.mat[1:13,14:26]
head(p.mat)


#setEPS()
#postscript("corrplot.eps")
pdf("corrplot.pdf")
#tiff("corrplot.tiff", height=6, width=6, units="in", res=600)
cex.before <- par("cex")
par(cex=0.7)
col <- colorRampPalette(c("#BB4444", "#EE9988", "#FFFFFF", "#77AADD", "#4477AA"))
corrplot(correlation, method="color", col=col(200),
type="full", order="original",
addCoef.col="black",
tl.col="black", tl.srt=45, tl.cex=1,
cl.cex=1/par("cex"),
p.mat=p.mat, sig.level=0.10, insig="blank",
addgrid.col="grey",
diag=TRUE
)
par(cex=cex.before)
dev.off()


######################################################
#############  Time Series Analysis  #################
######################################################

rm(list=ls())

timedata<-read.dta("time_series_data.dta")
attach(timedata)

#Create time series objects for each variable
t1nyt<-ts(nyt1, start=1994, frequency=1)
t2nyt<-ts(nyt2, start=1994, frequency=1)
t3nyt<-ts(nyt3, start=1994, frequency=1)
t4nyt<-ts(nyt4, start=1994, frequency=1)
t5nyt<-ts(nyt5, start=1994, frequency=1)
t6nyt<-ts(nyt6, start=1994, frequency=1)
t7nyt<-ts(nyt7, start=1994, frequency=1)
t8nyt<-ts(nyt8, start=1994, frequency=1)
t10nyt<-ts(nyt10, start=1994, frequency=1)
t12nyt<-ts(nyt12, start=1994, frequency=1)
t13nyt<-ts(nyt13, start=1994, frequency=1)
t14nyt<-ts(nyt14, start=1994, frequency=1)
t15nyt<-ts(nyt15, start=1994, frequency=1)
t1hearing<-ts(hearing1, start=1994, frequency=1)
t2hearing<-ts(hearing2, start=1994, frequency=1)
t3hearing<-ts(hearing3, start=1994, frequency=1)
t4hearing<-ts(hearing4, start=1994, frequency=1)
t5hearing<-ts(hearing5, start=1994, frequency=1)
t6hearing<-ts(hearing6, start=1994, frequency=1)
t7hearing<-ts(hearing7, start=1994, frequency=1)
t8hearing<-ts(hearing8, start=1994, frequency=1)
t10hearing<-ts(hearing10, start=1994, frequency=1)
t12hearing<-ts(hearing12, start=1994, frequency=1)
t13hearing<-ts(hearing13, start=1994, frequency=1)
t14hearing<-ts(hearing14, start=1994, frequency=1)
t15hearing<-ts(hearing15, start=1994, frequency=1)
t1panews<-ts(panews1, start=1994, frequency=1)
t2panews<-ts(panews2, start=1994, frequency=1)
t3panews<-ts(panews3, start=1994, frequency=1)
t4panews<-ts(panews4, start=1994, frequency=1)
t5panews<-ts(panews5, start=1994, frequency=1)
t6panews<-ts(panews6, start=1994, frequency=1)
t7panews<-ts(panews7, start=1994, frequency=1)
t8panews<-ts(panews8, start=1994, frequency=1)
t10panews<-ts(panews10, start=1994, frequency=1)
t12panews<-ts(panews12, start=1994, frequency=1)
t13panews<-ts(panews13, start=1994, frequency=1)
t14panews<-ts(panews14, start=1994, frequency=1)
t15panews<-ts(panews15, start=1994, frequency=1)


##### Plot of all time series data ####
#setEPS()
#postscript("ts_summary.eps")
pdf("ts_summary.pdf")
#tiff("ts_summary.tiff", height=6, width=6, units="in", res=600)
ts.plot(t1hearing, t2hearing, t3hearing, t4hearing, t5hearing, t6hearing, t7hearing, t8hearing, t10hearing, t12hearing, t13hearing, t14hearing, t15hearing, gpars=list(ylab="Percentage", xlab="Year", col=c("black"), lty=c(1), ylim=c(0,100), bty="n", axes=FALSE))
axis(1, at=c(1994:2005), labels=c(1994:2005))
axis(2, at=seq(0,100, by=10), labels=seq(0,100, by=10))
lines(t1nyt, col="red", lty=2)
lines(t2nyt, col="red", lty=2)
lines(t3nyt, col="red", lty=2)
lines(t4nyt, col="red", lty=2)
lines(t5nyt, col="red", lty=2)
lines(t6nyt, col="red", lty=2)
lines(t7nyt, col="red", lty=2)
lines(t8nyt, col="red", lty=2)
lines(t10nyt, col="red", lty=2)
lines(t12nyt, col="red", lty=2)
lines(t13nyt, col="red", lty=2)
lines(t14nyt, col="red", lty=2)
lines(t15nyt, col="red", lty=2)
lines(t1panews, col="blue", lty=5)
lines(t2panews, col="blue", lty=5)
lines(t3panews, col="blue", lty=5)
lines(t4panews, col="blue", lty=5)
lines(t5panews, col="blue", lty=5)
lines(t6panews, col="blue", lty=5)
lines(t7panews, col="blue", lty=5)
lines(t8panews, col="blue", lty=5)
lines(t10panews, col="blue", lty=5)
lines(t12panews, col="blue", lty=5)
lines(t13panews, col="blue", lty=5)
lines(t14panews, col="blue", lty=5)
lines(t15panews, col="blue", lty=5)
text(1994, 55, c("Commerce"), cex=0.8, col="red", pos=4)
legend(1994, 95, legend=c("Hearings", "NYT", "PA News"), lty=c(1,2,5), col=c("black", "red", "blue"), bty="n")
dev.off()

#############################################################################
############################# Figure 4 ######################################
#############################################################################


#setEPS()
#postscript("ts_1.eps")
pdf("ts_1.pdf", height=11, width=8.5)
#tiff("ts_1.tiff", height=11, width=8.5, units="in", res=600)

layout.matrix <- rbind(c(1,1,2,2), c(3,3,4,4), c(5,5,6,6)) #Create matrix for laying out plots
layout(layout.matrix) #Creates grid in plot window based on the layout matrix
#layout.show(6) #uncomment to show the seven panels in your plot window. Useful for re-arranging plots

par(mar=c(5,4,1,1))

#Economic
ts.plot(t1hearing, t1nyt, t1panews, gpars=list(ylab="", xlab="", col=c("black", "red", "blue"), lty=c(1,2,5), bty="n", axes=F, ylim=c(0,12), lwd=1.5))
title(main="", xlab="(a) Fiscal and Economic Issues", ylab="Percentage", cex=1.5, font.lab=2)
axis(1, at=c(1994:2005), labels=c(1994:2005), cex=1.5, font=2)
axis(2, at=seq(0,12, by=2), labels=seq(0,12, by=2), cex=1.5, font=2)
#minor.tick(nx=1, nx=1, tick.ratio=0.5)
legend("topleft", legend=c("House Hearings", "NYT", "PA News"), lty=c(1,2,5), col=c("black", "red", "blue"), bty="n", text.font=2, lwd=2)

#Civil Rights
ts.plot(t2hearing, t2nyt, t2panews, gpars=list(ylab="", xlab="", col=c("black", "red", "blue"), lty=c(1,2,5), bty="n", axes=F, lwd=1.5))
title(main="", xlab="(b) Civil Rights and Liberties", ylab="Percentage", cex=1.5, font.lab=2)
axis(1, at=c(1994:2005), labels=c(1994:2005), cex=1.5, font=2)
axis(2, at=seq(0,6, by=1), labels=seq(0,6, by=1), cex=1.5, font=2)
#minor.tick(nx=1, nx=1, tick.ratio=0.5)
#legend("topleft", legend=c("Hearings", "NYT", "PA News"), lty=c(1,2,5), col=c("black", "red", "blue"), bty="n")

#Health
ts.plot(t3hearing, t3nyt, t3panews, gpars=list(ylim=c(0,24), ylab="", xlab="", col=c("black", "red", "blue"), lty=c(1,2,5), bty="n", axes=F, lwd=1.5))
title(main="", xlab="(c) Health", ylab="Percentage", cex=1.5, font.lab=2)
axis(1, at=c(1994:2005), labels=c(1994:2005), cex=1.5, font=2)
axis(2, at=seq(0,24, by=2), labels=seq(0,24, by=2), cex=1.5, font=2)
#minor.tick(nx=1, nx=1, tick.ratio=0.5)
#legend("topleft", legend=c("Hearings", "NYT", "PA News"), lty=c(1,2,5), col=c("black", "red", "blue"), bty="n")

#Agriculture
ts.plot(t4hearing, t4nyt, t4panews, gpars=list(ylim=c(0,6), ylab="", xlab="", col=c("black", "red", "blue"), lty=c(1,2,5), bty="n", axes=F, lwd=1.5))
title(main="", xlab="(d) Agriculture", ylab="Percentage", cex=1.5, font.lab=2)
axis(1, at=c(1994:2005), labels=c(1994:2005), cex=1.5, font=2)
axis(2, at=seq(0,6, by=1), labels=seq(0,6, by=1), cex=1.5, font=2)
#minor.tick(nx=1, nx=1, tick.ratio=0.5)
#legend("topleft", legend=c("Hearings", "NYT", "PA News"), lty=c(1,2,5), col=c("black", "red", "blue"), bty="n")

#Labor
ts.plot(t5hearing, t5nyt, t5panews, gpars=list(ylim=c(0,8), ylab="", xlab="", col=c("black", "red", "blue"), lty=c(1,2,5), bty="n", axes=F, lwd=1.5))
title(main="", xlab="(e) Labor, Employment, Immigration", ylab="Percentage", cex=1.5, font.lab=2)
axis(1, at=c(1994:2005), labels=c(1994:2005), cex=1.5, font=2)
axis(2, at=seq(0,8, by=1), labels=seq(0,8, by=1), cex=1.5, font=2)
#minor.tick(nx=1, nx=1, tick.ratio=0.5)
#legend("topleft", legend=c("Hearings", "NYT", "PA News"), lty=c(1,2,5), col=c("black", "red", "blue"), bty="n")

#Education
ts.plot(t6hearing, t6nyt, t6panews, gpars=list(ylim=c(0,25), ylab="", xlab="", col=c("black", "red", "blue"), lty=c(1,2,5), bty="n", axes=F, lwd=1.5))
title(main="", xlab="(f) Education", ylab="Percentage", cex=1.5, font.lab=2)
axis(1, at=c(1994:2005), labels=c(1994:2005), cex=1.5, font=2)
axis(2, at=seq(0,24, by=2), labels=seq(0,24, by=2), cex=1.5, font=2)
#minor.tick(nx=1, nx=1, tick.ratio=0.5)
#legend("topleft", legend=c("Hearings", "NYT", "PA News"), lty=c(1,2,5), col=c("black", "red", "blue"), bty="n")

dev.off()


#setEPS()
#postscript("ts_2.eps")
#pdf("ts_2.pdf", height=11, width=8.5)
tiff("ts_2.tiff", height=11, width=8.5, units="in", res=600)
layout.matrix <- rbind(c(1,1,2,2), c(3,3,4,4), c(5,5,6,6), c(8,7,7,9)) #Create matrix for laying out plots
layout(layout.matrix) #Creates grid in plot window based on the layout matrix
#layout.show(9) #uncomment to show the seven panels in your plot window. Useful for re-arranging plots

par(mar=c(5,4,1,1))

#Environment
ts.plot(t7hearing, t7nyt, t7panews, gpars=list(ylim=c(0,12), ylab="", xlab="", col=c("black", "red", "blue"), lty=c(1,2,5), bty="n", axes=F, lwd=1.5))
title(main="", xlab="(g) Environment", ylab="Percentage", cex=1.5, font.lab=2)
axis(1, at=c(1994:2005), labels=c(1994:2005), cex=1.5, font=2)
axis(2, at=seq(0,12, by=2), labels=seq(0,12, by=2), cex=1.5, font=2)
#minor.tick(nx=1, nx=1, tick.ratio=0.5)
#legend("topleft", legend=c("Hearings", "NYT", "PA News"), lty=c(1,2,5), col=c("black", "red", "blue"), bty="n")

#Energy
ts.plot(t8hearing, t8nyt, t8panews, gpars=list(ylim=c(0,8), ylab="", xlab="", col=c("black", "red", "blue"), lty=c(1,2,5), bty="n", axes=F, lwd=1.5))
title(main="", xlab="(h) Energy", ylab="Percentage", cex=1.5, font.lab=2)
axis(1, at=c(1994:2005), labels=c(1994:2005), cex=1.5, font=2)
axis(2, at=seq(0,8, by=1), labels=seq(0,8, by=1), cex=1.5, font=2)
#minor.tick(nx=1, nx=1, tick.ratio=0.5)
#legend("topleft", legend=c("Hearings", "NYT", "PA News"), lty=c(1,2,5), col=c("black", "red", "blue"), bty="n")

#Transportation
ts.plot(t10hearing, t10nyt, t10panews, gpars=list(ylim=c(0,16), ylab="", xlab="", col=c("black", "red", "blue"), lty=c(1,2,5), bty="n", axes=F, lwd=1.5))
title(main="", xlab="(i) Transportation", ylab="Percentage", cex=1.5, font.lab=2)
axis(1, at=c(1994:2005), labels=c(1994:2005), cex=1.5, font=2)
axis(2, at=seq(0,16, by=4), labels=seq(0,16, by=4), cex=1.5, font=2)
#minor.tick(nx=1, nx=1, tick.ratio=0.5)
#legend("topleft", legend=c("Hearings", "NYT", "PA News"), lty=c(1,2,5), col=c("black", "red", "blue"), bty="n")

#Crime
ts.plot(t12hearing, t12nyt, t12panews, gpars=list(ylim=c(0,26), ylab="", xlab="", col=c("black", "red", "blue"), lty=c(1,2,5), bty="n", axes=F, lwd=1.5))
title(main="", xlab="(j) Law, Crime, and Family", ylab="Percentage", cex=1.5, font.lab=2)
axis(1, at=c(1994:2005), labels=c(1994:2005), cex=1.5, font=2)
axis(2, at=seq(0,26, by=4), labels=seq(0,26, by=4), cex=1.5, font=2)
#minor.tick(nx=1, nx=1, tick.ratio=0.5)
#legend("topleft", legend=c("Hearings", "NYT", "PA News"), lty=c(1,2,5), col=c("black", "red", "blue"), bty="n")

#Welfare
ts.plot(t13hearing, t13nyt, t13panews, gpars=list(ylim=c(0,10), ylab="", xlab="", col=c("black", "red", "blue"), lty=c(1,2,5), bty="n", axes=F, lwd=1.5))
title(main="", xlab="(k) Social Welfare", ylab="Percentage", cex=1.5, font.lab=2)
axis(1, at=c(1994:2005), labels=c(1994:2005), cex=1.5, font=2)
axis(2, at=seq(0,10, by=2), labels=seq(0,10, by=2), cex=1.5, font=2)
#minor.tick(nx=1, nx=1, tick.ratio=0.5)
#legend("topleft", legend=c("Hearings", "NYT", "PA News"), lty=c(1,2,5), col=c("black", "red", "blue"), bty="n")

#Housing
ts.plot(t14hearing, t14nyt, t14panews, gpars=list(ylim=c(0,16), ylab="", xlab="", col=c("black", "red", "blue"), lty=c(1,2,5), bty="n", axes=F, lwd=1.5))
title(main="", xlab="(l) Community Development and Housing", ylab="Percentage", cex=1.5, font.lab=2)
axis(1, at=c(1994:2005), labels=c(1994:2005), cex=1.5, font=2)
axis(2, at=seq(0,16, by=2), labels=seq(0,16, by=2), cex=1.5, font=2)
#minor.tick(nx=1, nx=1, tick.ratio=0.5)
#legend("topleft", legend=c("Hearings", "NYT", "PA News"), lty=c(1,2,5), col=c("black", "red", "blue"), bty="n")

#Commerce
ts.plot(t15hearing, t15nyt, t15panews, gpars=list(ylim=c(0,70), ylab="", xlab="", col=c("black", "red", "blue"), lty=c(1,2,5), bty="n", axes=F, lwd=1.5))
title(main="", xlab="(m) Banking, Finance, and Commerce", ylab="Percentage", cex=1.5, font.lab=2)
axis(1, at=c(1994:2005), labels=c(1994:2005), cex=1.5, font=2)
axis(2, at=seq(0,70, by=10), labels=seq(0,70, by=10), cex=1.5, font=2)
#minor.tick(nx=1, nx=1, tick.ratio=0.5)
#legend("topleft", legend=c("Hearings", "NYT", "PA News"), lty=c(1,2,5), col=c("black", "red", "blue"), bty="n")
dev.off()




